package math2;

/**
 * Gauss points data for Gauss-Legendre integration scheme on triangle
 * {-1<xi1,xi2<1, xi1+xi2<1}
 * 
 * @author hbui
 * 
 */
public class GaussTetrahedra {

	private static final double tet21 = (5 - Math.sqrt(5)) / 20;
	private static final double tet22 = (5 + 3 * Math.sqrt(5)) / 20;

	private static final double tet31 = (7 - Math.sqrt(15)) / 34;
	private static final double tet32 = (13 + 3 * Math.sqrt(15)) / 34;
	private static final double wtet31 = (2665 + 14 * Math.sqrt(15)) / 226800;
	private static final double tet33 = (7 + Math.sqrt(15)) / 34;
	private static final double tet34 = (13 - 3 * Math.sqrt(15)) / 34;
	private static final double wtet32 = (2665 - 14 * Math.sqrt(15)) / 226800;
	private static final double tet35 = (10 - 2 * Math.sqrt(15)) / 40;
	private static final double tet36 = (10 + 2 * Math.sqrt(15)) / 40;

	public static final double[][][] GAUSS_PTS = { //
	{ { 0.25, 0.25, 0.25, 1.00 / 6 } }, //
			{ { tet21, tet21, tet21, 1.00 / 24 }, //
					{ tet22, tet21, tet21, 1.00 / 24 }, //
					{ tet21, tet22, tet21, 1.00 / 24 }, //
					{ tet21, tet21, tet22, 1.00 / 24 } },
			// { { 0.25, 0.25, 0.25, -0.1333333333333333333333333333333 },
			// { 0.5, 0.1666666666666667, 0.1666666666666667, 0.075 },
			// { 0.1666666666666667, 0.5, 0.1666666666666667, 0.075 },
			// { 0.1666666666666667, 0.1666666666666667, 0.5, 0.075 },
			// { 0.1666666666666667, 0.1666666666666667, 0.1666666666666667,
			// 0.075 } }
			{ { 0.25, 0.25, 0.25, 16.0 / 810 }, //
					{ tet31, tet31, tet31, wtet31 }, //
					{ tet32, tet31, tet31, wtet31 }, //
					{ tet31, tet32, tet31, wtet31 }, //
					{ tet31, tet31, tet32, wtet31 }, //

					{ tet33, tet33, tet33, wtet32 }, //
					{ tet34, tet33, tet33, wtet32 }, //
					{ tet33, tet34, tet33, wtet32 }, //
					{ tet33, tet33, tet34, wtet32 }, //

					{ tet35, tet35, tet36, 20.0 / 2268 }, //
					{ tet35, tet36, tet35, 20.0 / 2268 }, //
					{ tet36, tet35, tet35, 20.0 / 2268 }, //
					{ tet36, tet36, tet35, 20.0 / 2268 }, //
					{ tet36, tet35, tet36, 20.0 / 2268 }, //
					{ tet35, tet36, tet36, 20.0 / 2268 } //
			} //
	};

	/**
	 * get number of Gauss points for a specific order
	 * 
	 * @param p
	 *            order
	 * @return number of Gauss points associated with order i
	 */
	public static int getNumGaussPoints(int p) {
		if (p > GAUSS_PTS.length) {
			throw new Error(
					"exceed maximum order of Gauss integration rule on tetrahedra. The maximum order is "
							+ p + " only");
		}
		return GAUSS_PTS[p - 1].length;
	}

	/**
	 * get xi1 coordinate of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return xi1
	 */
	public static double getGaussXi1(int p, int i) {
		return GAUSS_PTS[p - 1][i][0];
	}

	/**
	 * get xi2 coordinate of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return xi2
	 */
	public static double getGaussXi2(int p, int i) {
		return GAUSS_PTS[p - 1][i][1];
	}

	/**
	 * get xi2 coordinate of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return xi2
	 */
	public static double getGaussXi3(int p, int i) {
		return GAUSS_PTS[p - 1][i][2];
	}

	/**
	 * get weight of Gauss point j of order i
	 * 
	 * @param p
	 *            order of Gauss integration
	 * @param i
	 *            index of Gauss point
	 * @return weight
	 */
	public static double getGaussWeight(int p, int i) {
		return GAUSS_PTS[p - 1][i][3];
	}
}
